Vibration analysis method, program, storage medium

ABSTRACT

This vibration analysis method is a method for analyzing vibrations of a large-scale system with local strong nonlinearities, and includes a process (1) of applying the new type of complex modal analysis to an equation for a linear state variable to convert the equation to a real modal equation for lower-order modes, and correcting an effect of higher-order modes of the linear state variable from an equation for a nonlinear state variable and eliminating the modes, a process (2) of selecting secondary modes, which have a large effect on a solution of an original large-scale system, from the real modal equation for lower-order modes, and, in relation to secondary modes, which have a small effect, eliminating the modes thereof by incorporating the effect to the equation for nonlinear state variables as a correction term obtained from an approximate solution of the real modal equation for lower-order modes, and deriving the dimension reduced model, and a process (3) of calculating a frequency response by using the dimension reduced model.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority based on Japanese Patent ApplicationNo. 2020-161727 filed in Japan on Sep. 28, 2020, the content of which isincorporated herein.

TECHNICAL FIELD

The present invention mainly relates to a high-performance vibrationanalysis method for large-scale systems with local strongnonlinearities. More specifically, it relates to a dimensionalityreduction method using a new type of complex modal analysis.

BACKGROUND

A machine is composed of a large number of elements that move relativeto each other and has elements with strong nonlinearities such asbearings at base supporting portions or gears at joints between theelements. Moreover, in stress analysis and vibration analysis at designstage, by using a finite element method, and the like, it is modeled asa precise large-scale degree-of-freedom system close to an actualmachine in many cases. Therefore, almost all machines are large-scaledegree-of-freedom systems with local strong nonlinearities.

On the other hand, in recent years, there has been an increasing demandfor higher performance of machines. To realize this requirement, it isnecessary to implement a vibration analysis that fully considers aneffect of nonlinearity that is unavoidably included in the system at adesign stage. Moreover, to clarify the full picture of nonlinearvibrations that are complicated and diverse, it is not sufficient tosimply perform a numerical simulation for a specific system parameter,and the full content of frequency response including the stabilityanalysis needs to be clarified.

Until now, a Galerkin method based on linear eigenmode expansion (referto, for example, Patent Document 1 described below) has often been usedfor a vibration analysis of multi-degree-of-freedom nonlinear systemsbecause it has an algorithm with regularity and can be applied to a widerange of problems.

SUMMARY

However, the Galerkin method has significant disadvantages such as anexplosive increase in calculation time and memory size when an analysiswith high accuracy is attempted by increasing the number of modes to beused and the disappearance of a mode separation effect with nonlinearityby dispersing throughout a system after a conversion to modalcoordinates even if a model has nonlinearity only locally.

In order to overcome such a situation, research has been vigorouslyconducted to perform an analysis with high accuracy using a dimensionreduced model with a small number of modes but have not yet reached afundamental solution. The main problem is that accuracy of the analysis(especially an accuracy of stability analysis) deterioratessignificantly due to a coupling between modes specific to a nonlinearsystem if modes to be used are not appropriately selected. In thismanner, it is not an exaggeration to say that there are no practicalvibration analysis methods for a large-scale nonlinear system, and it isan urgent problem in numerical calculation related to mechanical designto develop a high-performance vibration analysis method that can dealwith this problem.

The present invention has been made with aims to fundamentally solve theproblems described above and provides a method of rationallyconstructing the dimension reduced model with high accuracy in which theeffect of nonlinearity is appropriately reflected for a large-scalesystem with local strong nonlinearity.

A proposed method is configured mainly from the following two processes.(1) A process of applying the new type of complex modal analysis to theequations for linear state variables [Equation (1) and Equation (4) tobe described below] to convert the equations to the real modal equationsfor lower-order modes, and correcting and eliminating an effect ofhigher-order modes of the linear state variables from the equation fornonlinear state variables [Equation (2) to be described below], and (2)a process of selecting modes (dominant modes), which have a large effecton a solution of an original large-scale system, from the real modalequation for lower-order modes, and, in relation to modes (secondarymodes), which have a small effect, eliminating secondary modes byincorporating the effect to the equation for the nonlinear statevariables [Equation (2)] as correction terms obtained from anapproximate and deriving the dimension reduced model, and a process ofselecting dominant modes from lower-order modes (as rationally aspossible) in a calculation procedure of a frequency response

By using this dimension reduced model, a vibration analysis includingstability analysis of a large-scale nonlinear system, which has not beenpossible, can be performed at high speed and with high accuracy.Moreover, the proposed method has very high versatility and can begenerally applied to large-scale vibration systems with local strongnonlinearities.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram which shows an overview of the present analysistechnique.

FIG. 2 is a flowchart of the present analysis technique.

FIG. 3 is a schematic diagram which shows an analytical model of anapplication example of a dimensionality reduction method.

FIG. 4 is a frequency response diagram which shows first-order tofifth-order primary resonance regions in the full model.

FIG. 5 is a frequency response diagram which shows the first-order andsecond-order primary resonance regions in the full model.

FIG. 6 is a frequency response diagram which shows the first-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=59, n_(a)=5).

FIG. 7 is a frequency response diagram which shows the first-order andsecond-order primary resonance regions in the dimension reduced model(n_(l)=59, n_(a)=5).

FIG. 8 is a frequency response diagram which shows the first-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=59, n_(a)=15).

FIG. 9 is a frequency response diagram of the first-order andsecond-order primary resonance regions in the dimension reduced model(n_(l)=59, n₈=15).

FIG. 10 is a frequency response diagram which shows the first-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=5, n_(a)=5).

FIG. 11 is a frequency response diagram which shows the first-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=10, n_(a)=5).

FIG. 12 is a diagram which shows an effect of the lower-order modesnumber n_(l) on the first-order natural frequency.

FIG. 13 is a diagram which shows an effect of the lower-order modesnumber n_(l) on the second-order natural frequency.

FIG. 14 is a diagram which shows an effect of the lower-order modesnumber n_(l) on the third-order natural frequency.

FIG. 15 is a diagram which shows an effect of the lower-order modesnumber n_(l) on the fourth-order natural frequency.

FIG. 16 is a diagram which shows an effect of the lower-order modesnumber n_(l) on the fifth-order natural frequency.

FIG. 17 is a frequency response diagram of the second-order primaryresonance region in the full model.

FIG. 18 is a frequency response diagram of the second-order primaryresonance region in the dimension reduced model (n_(l)=20, δ=5×10⁻⁴).

FIG. 19 is a frequency response diagram of the second-order primaryresonance region in the dimension reduced model (n_(l)=20, δ=1×10⁻⁴).

FIG. 20 is a frequency response diagram of the second-order primaryresonance region in the dimension reduced model (n_(l)=20, δ=1×10⁻⁵).

FIG. 21 is a frequency response diagram of the first-order tosecond-order primary resonance regions in the dimension reduced model(n_(l)=15, δ=7×10⁻⁴).

FIG. 22 is a frequency response diagram of the first-order tosecond-order primary resonance regions in the dimension reduced model(n_(l)=15, δ=1×10⁻⁴).

FIG. 23 is a frequency response diagram of the first-order tosecond-order primary resonance regions in the dimension reduced model(n_(l)=15, δ=2×10⁻⁵).

FIG. 24 is a frequency response diagram of the first-order tosecond-order primary resonance regions in the dimension reduced model(n_(l)=599, n_(a)=20).

FIG. 25 is a frequency response diagram of the fourth-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=30, δ=1×10⁻⁴).

FIG. 26 is a frequency response diagram of the fourth-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=30, δ=5×10⁻⁴).

FIG. 27 is the frequency response diagram of the fourth-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=30, δ=2×10⁻⁴).

FIG. 28 is a frequency response diagram of the fourth-order tofifth-order primary resonance regions in the dimension reduced model(n_(l)=599, n_(a)=20).

FIG. 29 is a block diagram which shows an example of a device thatimplements the present analysis technique.

DETAILED DESCRIPTION

Specific implementation items in the two processes described above willbe described below.

It should be noted that, in the present specification, expressions ofterms and mathematical formulas are restricted in the sentences otherthan independent mathematical formulas that are captured as the images.For this reason, for example, a term captured as an image in the presentspecification and expressed as a term (0-1) is expressed as “^(C)ψ{tildeover ( )}₁ ^(T)” in the sentence of the present specification. A termthat is incorporated as an image in the present specification andexpressed as a term (0-2) is expressed as “y{dot over ( )}” in the textof the present specification. Note that “y{dot over ( )}” means a termobtained by first-order differentiation of y with respect to time.

^(C){tilde over (ψ)}₁ ^(T)  (0-1)

{dot over (y)}  (0-2)

A method proposed in the present invention analyzes a large-scalevibration system with N degrees of freedom and local strongnonlinearities. The fundamental equations can be generally expressed asfollows.

$\begin{matrix}{{{{M_{11}{\overset{¨}{y}}_{1}} + {C_{11}{\overset{.}{y}}_{1}} + {K_{11}y_{1}} + {M_{12}{\overset{¨}{y}}_{2}} + {C_{12}{\overset{.}{y}}_{2}} + {K_{12}y_{2}}} = f_{1}},} & (1)\end{matrix}$ $\begin{matrix}{{{M_{22}{\overset{¨}{y}}_{2}} + {C_{22}{\overset{.}{y}}_{2}} + {K_{22}y_{2}} + {M_{21}{\overset{¨}{y}}_{1}} + {C_{21}{\overset{.}{y}}_{1}} + {K_{21}y_{1}} + {n_{2}\left( {y_{2},{\overset{.}{y}}_{2}} \right)}} = f_{2}} & (2)\end{matrix}$ Where, $\begin{matrix}\left. \begin{matrix}{y_{1},{f_{1} = {{f_{1}\left( {\omega t} \right)}:\left( {n \times 1} \right)}},y_{2},{n_{2} = {n_{2}\left( {y_{2},{\overset{.}{y}}_{2}} \right)}},{f_{2} = {{f_{2}\left( {\omega t} \right)}:\left( {m \times 1} \right)}},{\text{?} = {d/{dt}}}} \\{M_{11},C_{11},{K_{11}:\left( {n \times n} \right)},M_{12},C_{12},{K_{12}:\left( {n \times m} \right)}} \\{M_{22},C_{22},{K_{22}:\left( {m \times m} \right)},M_{21},C_{21},{K_{21}:\left( {m \times n} \right)},{{n + m} = N},{m\text{?}n}}\end{matrix} \right\} & (3)\end{matrix}$ ?indicates text missing or illegible when filed

y₁ is an n-dimensional vector that summarizes the physical coordinatesof state variables on which the nonlinear forces do not act, and y₂ isan m-dimensional vector that summarizes the physical coordinates ofstate variables on which the nonlinear forces act. In the followingdescription, y₁ will be referred to as a linear state variable and y₂will be referred to as a nonlinear state variable. For convenience,Equation (1) will be referred to as a linear state equation and Equation(2) will be referred to as a nonlinear state equation. n₂ is set to astrong nonlinear function with respect to y₂ and y{dot over ( )}2, andf₁ and f₂ are set to the harmonic forced external forces of a period2π/ω. M is a coefficient matrix (a mass matrix) for mass elements. C isa coefficient matrix (a damping matrix) for damping elements. K is acoefficient matrix (a stiffness matrix) for linear spring elements.Positive definiteness for M₁₁, C₁₁ and K₁₁ are not assumed to ensureversatility of the proposed method. The fact that m<<<n indicates thatthe nonlinearity exists only locally. In this manner, the proposedmethod is applicable to almost large-scale vibration systems with localstrong nonlinearities.

In formulating a dimensionality reduction method, a following system inwhich the coefficient matrices of Equation (1) are replaced with thetransposed matrices and the external force set to zero is treated at thesame time.

M ₁₁ ^(T) ÿ ₁ *+C ₁₁ ^(T) {dot over (y)} ₁ *+K ₁₁ ^(T) y ₁ *+M ₂₁ ^(T) ÿ₂ *C ₂₁ ^(T) {dot over (y)} ₂ *+K ₂₁ ^(T) y ₂*=0  (4)

Where, an upper right subscript “T” is the transposition symbol. y₂ inEquation (4) is a solution that satisfies Equations (1) and (2), butsince y₁ is generally different from a solution that satisfies Equations(1) and (2), it is expressed as y*₁ in Equation (4). In the followingdescription, Equation (4) is referred to as a dual system of Equation(1).

To formulate the dimensionality reduction method based on the new typeof complex modal analysis, Equations (1) and (4) are transformed intothe first-order differential equations as follows. Note that In denotesan nth order unit matrix.

$\begin{matrix}{{{B_{11}{\overset{.}{x}}_{1}} + {B_{12}{\overset{.}{x}}_{2}}} = {{A_{11}x_{1}} + {A_{12}x_{2}} - {E_{1}g_{1}}}} & (5)\end{matrix}$ $\begin{matrix}{{{B_{11}^{T}{\overset{.}{x}}_{1}^{*}} + {B_{21}^{T}{\overset{.}{x}}_{2}}} = {{A_{11}^{T}x_{1}^{*}} + {A_{12}^{*}x_{2}}}} & (6)\end{matrix}$ Where, $\begin{matrix}\left. \begin{matrix}{{A_{11} = \begin{bmatrix}0 & K_{11} \\K_{11} & C_{11}\end{bmatrix}},{A_{12} = \begin{bmatrix}0 & 0 \\K_{12} & C_{12}\end{bmatrix}},{B_{11} = \begin{bmatrix}K_{11} & 0 \\0 & {- M_{11}}\end{bmatrix}},{B_{12} = \begin{bmatrix}0 & 0 \\0 & {- M_{12}}\end{bmatrix}}} \\{{A_{12}^{*} = \begin{bmatrix}0 & 0 \\K_{21}^{T} & C_{21}^{T}\end{bmatrix}},{B_{21} = \begin{bmatrix}0 & 0 \\0 & {- M_{21}}\end{bmatrix}},{E_{1} = \begin{bmatrix}0 & 0 \\\text{?} & 0\end{bmatrix}}} \\{{x_{1} = \begin{bmatrix}y_{1} \\{\overset{.}{y}}_{1}\end{bmatrix}},{x_{2} = \begin{bmatrix}y_{2} \\{\overset{.}{y}}_{2}\end{bmatrix}},{x_{1}^{*} = \begin{bmatrix}y_{1}^{*} \\{\overset{.}{y}}_{1}^{*}\end{bmatrix}},{g_{1} = \begin{bmatrix}f_{1} \\{\overset{.}{f}}_{1}\end{bmatrix}}}\end{matrix} \right\} & (7)\end{matrix}$ ?indicates text missing or illegible when filed

In the dimensionality reduction method, first, the nonlinear statevariables are fixed by setting x·2=0 and x₂=0 in Equations (5) and (6),and a following general eigenvalue problem is derived from a freevibration system with g₁=0.

[A ₁₁ −λB ₁₁ ]x ₁=0  (8)

[A ₁₁ ^(T) −λB ₁₁ ^(T) ]x ₁*=0  (9)

The eigenvalues λ obtained from Equations (8) and (9) are consistentwith each other and generally become complex numbers. In addition,corresponding eigenvectors X₁ and X*₁ are also complex vectors. X₁ is aright complex eigenvector and X*₁ is a left complex eigenvector forEquation (8), and X*₁ is a right complex eigenvector and X₁ is a leftcomplex eigenvector for Equation (9). In terms of vibration engineering,an imaginary part of the eigenvalue λ corresponds to a natural angularfrequency, and X₁ and X*₁ correspond to left and right complexconstraint modes.

As a first step, the complex eigenvalues and the left and right complexeigenvectors (the complex constraint modes) are obtained from Equations(8) and (9). However, when a dimension 2n of the general eigenvalueproblem is large, it is difficult to obtain all eigenvalues andeigenvectors (a solution of a large-scale general eigenvalue problemremains an unsolved conundrum in a field of numerical calculation and isvirtually almost impossible). However, it is still possible to obtaineigenvectors corresponding to a relatively small number of eigenvaluesin order from a lower order (in ascending order of an absolute value).

On the other hand, it is known that an effect on vibrationcharacteristics generally decreases as a mode order increases in bothlinear and nonlinear systems.

Therefore, in this dimensionality reduction method, only relativelylower-order eigenvalues and eigenvectors from the first order to n_(l)^(th) (n_(l)<<n) order, which may have a large effect on an accuracy ofa solution, are obtained from the general eigenvalue problem, an effectof higher-order modes of an (n₁+1)^(th) order or higher is approximatedusing lower-order modes up to the n_(l) ^(th) order. A value of n_(l)can be set based on a maximum order of a natural frequency (an imaginarypart of an eigenvalue) present within a frequency range about severaltimes a maximum frequency to be analyzed.

The relational equations required for this procedure are as follows. Inthe equations, subscripts “l” and “h” indicate variables or physicalquantities associated with lower-order modes and higher-order modes,respectively. A physical quantity with an upper left subscript “C”indicates that it is a complex number, and similarly, “R” and “l” inprinciple represent the real part and the imaginary part thereof,respectively.

First, left and right eigenvectors corresponding to eigenvalues from thefirst order to the n_(l) ^(th) (n₁<<n) order in Equation (8) areobtained. These eigenvalues are assumed to have negative complex numbersfor all real parts, and a p^(th) order (p=1, . . . , n_(l)) complexconjugate eigenvalue is expressed as follows.

$\begin{matrix}{\left. \begin{matrix}{\text{?},{\text{?} = {{\text{?} \pm \text{?}} = {{{- \zeta_{p,1}}\omega_{p,1}} \pm {i\sqrt{1 - \zeta_{p,1}^{2}}\omega_{p,1}}}}},{i = \sqrt{- 1}}} \\{{\text{?} = {{- \zeta_{p,1}}\omega_{p,1}}},{\text{?} = {\sqrt{1 - \zeta_{p,1}^{2}}\omega_{p,1}}},{\omega_{p,1} > 0},{0 \leq \zeta_{p,1} < 0}}\end{matrix} \right\}\left( {{p = 1},\text{?}} \right)} & (10)\end{matrix}$ ?indicates text missing or illegible when filed

Where, ω_(p,1) and ζ_(p,1) represent an undamped natural angularfrequency and a damping ratio of a p^(th) order mode, respectively.Furthermore, these n_(l) pairs of eigenvalues (2n_(l) eigenvalues intotal) are put together and displayed in a matrix as follows.

$\begin{matrix}\left. \begin{matrix}{\text{?},{\text{?} = {\text{?} \pm \text{?}}},{\text{?} = {{- Z_{\mu}}\Omega_{\mu}}},{\text{?} = {\sqrt{\text{?} - Z_{\mu}^{2}}\Omega_{\mu}:\left( {n_{i} \times n_{i}} \right)}}} \\{{\text{?} = {{diag}\left\lbrack {\text{?}\text{?}\cdots\text{?}} \right\rbrack}},{\text{?} = {{{diag}\left\lbrack {\text{?}\text{?}\cdots\text{?}} \right\rbrack}:\left( {n_{i} \times n_{i}} \right)}}} \\{{\Omega_{\mu} = {{{diag}\left\lbrack {\omega_{1,1}\omega_{2,1}\cdots\text{?}} \right\rbrack}:\left( {n_{i} \times n_{i}} \right)}},{\omega_{p,1} = \sqrt{\text{?} + \text{?}}}} \\{{Z_{\mu} = {{{diag}\left\lbrack {\zeta_{1,1}\zeta_{2,1}\cdots\text{?}} \right\rbrack}:\left( {n_{i} \times n_{i}} \right)}},{\zeta_{p,1} = \frac{\text{?}}{\omega_{p,1}}}}\end{matrix} \right\} & (11)\end{matrix}$ ?indicates text missing or illegible when filed

Where, diag[ ] represents a diagonal matrix.

A right complex constraint modal matrix in which n_(l) sets (2n_(l)pieces) of the right complex eigenvectors of Equation (8) correspondingto this complex eigenvalue are put together is set as ^(C)Φ{circumflexover ( )}_(1l). Similarly, a left complex constraint modal matrix inwhich n_(l) sets (2n_(l) pieces) of left complex eigenvectors are puttogether is set as ^(C)Ψ{circumflex over ( )}_(1l). It is assumed that^(C)Φ{circumflex over ( )}_(1l) and ^(C)Ψ{circumflex over ( )}_(1l) arenormalized to satisfy the following equation.

$\begin{matrix}\left. \begin{matrix}{\text{?} = {{\begin{bmatrix}{\text{?} - \text{?}} & 0 \\0 & {\text{?} - \text{?}}\end{bmatrix}\begin{bmatrix}\text{?} & 0 \\0 & \text{?}\end{bmatrix}}:\left( {2\text{?} \times 2\text{?}} \right)}} \\{\text{?} = {\begin{bmatrix}{\text{?} - \text{?}} & 0 \\0 & {\text{?} - \text{?}}\end{bmatrix}:\left( {2\text{?} \times 2\text{?}} \right.}}\end{matrix} \right\} & (12)\end{matrix}$ ?indicates text missing or illegible when filed

At this time, ^(C)Φ{circumflex over ( )}_(1l) and ^(C)Ψ{circumflex over( )}_(1l) are expressed as follows. Note that ^(R)φ_(p,1) and^(l)φ_(p,1) are a real part and an imaginary part of the right complexeigenvector of the p^(th) order mode, respectively. In addition,^(R)ψ_(p,1) and ^(l)ψ_(p,1) are a real part and an imaginary part of theleft complex eigenvector of the p^(th) order mode, respectively.

$\begin{matrix}\left. \begin{matrix}{{\text{?} = \begin{bmatrix}\text{?} & \text{?} \\\text{?} & \text{?}\end{bmatrix}},{\text{?} = {\begin{bmatrix}\text{?} & \text{?} \\\text{?} & \text{?}\end{bmatrix}:\left( {2n \times 2\text{?}} \right)}}} \\{\text{?},{\text{?} = {\text{?} \pm \text{?}}},\text{?},{\text{?} = {\text{?} \pm {\text{?}:\left( {n \times \text{?}} \right)}}}} \\{{\text{?} = \left\lbrack {\text{?}\text{?}\cdots\text{?}} \right\rbrack},{\text{?} = {\left\lbrack {\text{?}\text{?}\cdots\text{?}} \right\rbrack:\left( {n \times \text{?}} \right)}}} \\{{\text{?} = \left\lbrack {\text{?}\text{?}\cdots\text{?}} \right\rbrack},{\text{?} = {\left\lbrack {\text{?}\text{?}\cdots\text{?}} \right\rbrack:\left( {n \times \text{?}} \right)}}}\end{matrix} \right\} & (13)\end{matrix}$ ?indicates text missing or illegible when filed

In order to speed up the calculation, a right real constraint modalmatrix ^(R)Φ{circumflex over ( )}_(1l) and a left real constraint modalmatrix ^(R)ψ{circumflex over ( )}_(1l) are derived from the right andleft complex constraint modal matrices ^(C)Φ{circumflex over ( )}_(1l)and ^(C)Ψ{circumflex over ( )}_(1l) for the lower-order modes asfollows.

$\begin{matrix}\left. \begin{matrix}{\text{?} = {{\text{?}\begin{bmatrix}\text{?} & \text{?} \\\text{?} & \text{?}\end{bmatrix}}^{- 1} = {\begin{bmatrix}\text{?} & \text{?} \\\text{?} & \text{?}\end{bmatrix}:\left( {2n \times 2\text{?}} \right)}}} \\{\text{?} = {{\text{?}\begin{bmatrix}\text{?} & \text{?} \\\text{?} & \text{?}\end{bmatrix}}^{- 1} = {\begin{bmatrix}\text{?} & \text{?} \\\text{?} & \text{?}\end{bmatrix}:\left( {2n \times 2\text{?}} \right)}}} \\{\left. \begin{matrix}{{\text{?} = {\text{?} - \text{?}}},{\text{?} = {\text{?} + \text{?}}}} \\{{\text{?} = {\text{?} - \text{?}}},{\text{?} = {\text{?} + \text{?}}}} \\{{\text{?} = \text{?}},{\text{?} = \text{?}},{\text{?} = \text{?}},{\text{?} = \text{?}}}\end{matrix} \right\}:\left( {n \times \text{?}} \right)}\end{matrix} \right\} & (14)\end{matrix}$ ?indicates text missing or illegible when filed

Physical coordinates x₁ and x*₁ of the linear state variables can berepresented by a sum of x_(1l) and x*_(1l) based on the lower-ordermodes and x_(1h) and X*_(1h) based on the higher-order modes as follows.Note that y_(1l) and y*_(1l) are the linear state variables for thelower-order modes, respectively. y_(1h) and y*_(1h), are the linearstate variables for the higher-order modes, respectively.

$\begin{matrix}{{x_{1} = {x_{1l} + x_{1h}}},{x_{1l} = {\begin{bmatrix}y_{1l} \\{\overset{.}{y}}_{1l}\end{bmatrix}:\left( {2n \times 1} \right)}},{x_{1h} = {\begin{bmatrix}y_{1h} \\{\overset{.}{y}}_{1h}\end{bmatrix}:\left( {2n \times 1} \right)}}} & (15)\end{matrix}$ $\begin{matrix}{{x_{1}^{*} = {x_{1l}^{*} + x_{1h}^{*}}},{x_{1l}^{*} = {\begin{bmatrix}y_{1l}^{*} \\{\overset{.}{y}}_{1l}^{*}\end{bmatrix}:\left( {2n \times 1} \right)}},{x_{1h}^{*} = {\begin{bmatrix}y_{1h}^{*} \\{\overset{.}{y}}_{1h}^{*}\end{bmatrix}:\left( {2n \times 1} \right)}}} & (16)\end{matrix}$

Among these, for x_(1l) and x*_(1l), real modal coordinates η_(1l) andη*_(1l) are introduced by the following equations via ^(C)Φ{circumflexover ( )}_(1l) and ^(C)Ψ{circumflex over ( )}_(1l) respectively. Notethat ξ_(1l) is the real modal coordinates corresponding for the physicalcoordinates of the lower-order modes.

$\begin{matrix}\left. \begin{matrix}{{x_{1l} =^{R}{{\overset{.}{\Phi}}_{1l}\left( {\eta_{1l} + {T_{12l}x_{2}} + {T_{11l}g_{1}}} \right)}},{\eta_{1l} = {\begin{bmatrix}\xi_{1l} \\{\overset{.}{\xi}}_{1l}\end{bmatrix}:\left( {2n_{l} \times 1} \right)}}} \\\begin{matrix}{{T_{12l} = {\begin{bmatrix}{{- {\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}}C_{12}} & {{- {\,^{J}\Psi_{1l}^{T}}}M_{12}} \\{{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}K_{12}} & 0\end{bmatrix}:\left( {2n_{l} \times 2m} \right)}},} \\{T_{11l} = {\begin{bmatrix}0 & 0 \\{- {\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}} & 0\end{bmatrix}:\left( {2n_{l} \times 2n} \right)}}\end{matrix}\end{matrix} \right\} & (17)\end{matrix}$ $\begin{matrix}\left. \begin{matrix}{{x_{1l}^{*} = {{\,^{R}{\overset{\sim}{\Psi}}_{1l}}\left( {\eta_{1l}^{*} + {T_{12l}^{*}x_{2}}} \right)}},{\eta_{1l}^{*} = {\begin{bmatrix}\xi_{1l}^{*} \\{\overset{.}{\xi}}_{1l}^{*}\end{bmatrix}:\left( {2n_{l} \times 1} \right)}}} \\{T_{12l}^{*} = {\begin{bmatrix}{{- {\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}}C_{21}^{T}} & {{- {\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}}M_{21}^{T}} \\{{\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}K_{21}^{T}} & 0\end{bmatrix}:\left( {2n_{l} \times 2m} \right)}}\end{matrix} \right\} & (18)\end{matrix}$

After substituting Equations (15) and (17) into Equation (5),multiplying it by ^(R)Ψ{circumflex over ( )}_(1l) ^(T) from the left andrearranging it, a real modal equation for the lower-order modes isderived as follows.

$\begin{matrix}\left. \begin{matrix}{{{\overset{¨}{\xi}}_{1l} + {2Z_{1l}\Omega_{1l}{\overset{.}{\xi}}_{1l}} + {\Omega_{1l}^{2}\xi_{1l}} + {R_{12l}{\overset{¨}{y}}_{2}} + {Q_{12l}{\overset{.}{y}}_{2}} + {P_{12l}y_{2}}} = {{{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}{\overset{.}{f}}_{1}} + {{\overset{\sim}{\Psi}}_{11l}^{T}f_{1}}}} \\\begin{matrix}{{P_{12l} = {{{- \Omega_{1l}^{2}}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}C_{12}} + {{\overset{\sim}{\Psi}}_{11l}^{T}K_{12}}}},} \\{Q_{12l} = {{{- \Omega_{1l}^{2}}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}M_{12}} + {{\overset{\sim}{\Psi}}_{22l}^{T}C_{12}} + {{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}K_{12}}}}\end{matrix} \\{{R_{12l} = {{\overset{\sim}{\Psi}}_{22l}^{T}M_{12}}},P_{12l},Q_{12l},{R_{12l}:\left( {n_{l} \times m} \right)}}\end{matrix} \right\} & (19)\end{matrix}$

For the dual system, substituting Equations (16) and (18) into Equation(6), multiplying it by ^(R)Φ{circumflex over ( )}_(1l) ^(T) from theleft and rearranging it, a real modal equation for the lower-order modesis derived as follows.

$\begin{matrix}\left. \begin{matrix}{{{\overset{¨}{\xi}}_{1l}^{*} + {2Z_{1l}\Omega_{1l}{\overset{.}{\xi}}_{1l}^{*}} + {\Omega_{1l}^{2}\xi_{1}^{*}} + {R_{12l}^{*}{\overset{¨}{y}}_{2}} + {Q_{12l}^{*}{\overset{.}{y}}_{2}} + {P_{12l}^{*}y_{2}}} = 0} \\\begin{matrix}{{P_{12l}^{*} = {{{- \Omega_{1l}^{2l}}{\overset{\sim}{\Phi}}_{1l}^{T}C_{21}^{T}} + {{\overset{\sim}{\Phi}}_{11l}^{T}K_{21}^{T}}}},} \\{Q_{12l}^{*} = {{{- \Omega_{1l}^{2l}}{\overset{\sim}{\Phi}}_{1l}^{T}M_{21}^{T}} + {{\overset{\sim}{\Phi}}_{22l}^{T}C_{21}^{T}} + {{\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}K_{21}^{T}}}}\end{matrix} \\{{R_{12l}^{*} = {{\overset{\sim}{\Phi}}_{22l}^{T}M_{21}^{T}}},P_{12l}^{*},Q_{12l}^{*},{R_{12l}^{*}:\left( {n_{l} \times m} \right)}}\end{matrix} \right\} & (20)\end{matrix}$

In this manner, even if the coefficient matrices do not satisfy positivedefiniteness, the mass, damping, and stiffness matrices can bediagonalized simultaneously, and the modal equations can be derived inthe form of real-type second-order differential equation. This is themost important feature of the new type of complex modal analysis.

The linear state variable y_(1l) in the nonlinear state equation[Equation (2)] is represented by a sum of y_(1l) based on thelower-order modes and y_(1h) based on the higher-order modes as shown inEquation (15), and y_(1l) is converted into the real modal coordinatescorresponding to the lower-order modes introduced in Equation (17).Furthermore, if the effect of y_(1h) based on the higher-order modes iscorrected and eliminated, the following equation is obtained.

$\begin{matrix}\left. \begin{matrix}{{\left( {M_{22} + {\overset{\sim}{M}}_{22{li}}} \right){\overset{¨}{y}}_{2}} + {\left( {C_{22} + {\overset{\sim}{C}}_{22{li}}} \right){\overset{.}{y}}_{2}} + {\left( {K_{22} + {\overset{\sim}{K}}_{22h}} \right)y_{2}} +} \\{{{R_{12l}^{*T}{\overset{¨}{\xi}}_{1l}} + {Q_{12l}^{*T}{\overset{.}{\xi}}_{1l}} + {P_{12l}^{*T}\xi_{1l}} + n_{2}} = {f_{2} + {q_{21h}{\overset{.}{f}}_{1}} + {p_{21h}f_{1}}}}\end{matrix} \right\} & (21)\end{matrix}$

Where, M{tilde over ( )}_(22h), C{tilde over ( )}_(22h), K{tilde over( )}_(22h), p_(21h), and q_(21h) represent correction terms of thehigher-order modes. These are all obtained from the lower-order modes asfollows.

$\begin{matrix}\left. \begin{matrix}\begin{matrix}{{\overset{\sim}{M}}_{22h} = {{M_{21}\Gamma_{1h}C_{12}} - {M_{21}\Gamma_{2h}K_{12}} + {C_{21}\Gamma_{1h}M_{12}} - {C_{21}\Gamma_{2h}C_{12}} -}} \\{{C_{21}\Gamma_{3h}K_{12}} - {K_{21}\Gamma_{2h}M_{12}} - {K_{21}\Gamma_{3h}C_{12}} - {K_{21}\Gamma_{4h}K_{12}}}\end{matrix} \\\begin{matrix}{{\overset{\sim}{C}}_{22h} = {{C_{21}\Gamma_{1h}C_{12}} - {K_{21}\Gamma_{2h}C_{12}} - {C_{21}\Gamma_{2h}K_{12}} -}} \\{{K_{21}\Gamma_{3h}K_{12}} + {M_{21}\Gamma_{1h}K_{12}}}\end{matrix} \\{{\overset{\sim}{K}}_{22h} = {{C_{21}\Gamma_{1h}K_{12}} - {K_{21}\Gamma_{2h}K_{12}}}} \\{p_{21h} = {{C_{21}\Gamma_{1h}} - {K_{21}\Gamma_{2h}}}} \\{q_{21h} = {{M_{21}\Gamma_{1h}} - {C_{21}\Gamma_{2h}} - {K_{21}\Gamma_{3h}}}} \\{\Gamma_{1h} = {{{{- {\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}}{\overset{\sim}{\Psi}}_{22l}^{T}} - {{\overset{\sim}{\Phi}}_{11l}{\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}}} = {{{- {\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}}{\overset{\sim}{\Psi}}_{11l}^{T}} - {{\overset{\sim}{\Phi}}_{22l}{\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}}}}} \\{\Gamma_{2h} = {K_{11}^{- 1} - {{\overset{\sim}{\Phi}}_{11l}\Omega_{1l}^{- 2}{\overset{\sim}{\Psi}}_{11l}^{T}} + {{\,^{l}{\overset{\sim}{\Phi}}_{1l}}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}}}} \\{\Gamma_{3h} = {{{- K_{11}^{- 1}}C_{11}K_{11}^{- 1}} - {2{\overset{\sim}{\Phi}}_{11l}^{R}\Lambda_{1l}\Omega_{1l}^{- l}{\overset{\sim}{\Psi}}_{11l}^{T}} - {{\overset{\sim}{\Phi}}_{11l}\Omega_{1l}^{{- 2}l}{\overset{\sim}{\Psi}}_{1l}^{T}} - {{\,^{l}{\overset{\sim}{\Phi}}_{1l}}\Omega_{1l}^{- 2}{\overset{\sim}{\Psi}}_{11l}^{T}}}} \\\begin{matrix}{\Gamma_{4h} = {{K_{11}^{- 1}C_{11}K_{11}^{- 1}C_{11}K_{11}^{- 1}} - {K_{11}^{- 1}M_{11}K_{11}^{- 1}} - {{\overset{\sim}{\Phi}}_{11l}\Omega_{1l}^{- 4}\left\{ \left( {{4\Omega_{1l}^{{- 2}R}\Lambda_{1l}^{2}} -} \right. \right.}}} \\{\left. {{\left. I_{u_{i}} \right){\overset{\sim}{\Psi}}_{11l}^{T}} + {2{\,^{R}\Lambda_{11}^{l}}{\overset{\sim}{\Psi}}_{11}^{T}}} \right\} - {{\,^{l}{\overset{\sim}{\Phi}}_{1l}}{\Omega_{1l}^{- 2}\left( {{2\Omega_{1l}^{- 2}{\,^{R}\Lambda_{1l}^{T}}{\overset{\sim}{\Psi}}_{11l}^{T}} + {\,^{l}{\overset{\sim}{\Phi}}_{1l}^{T}}} \right)}}}\end{matrix}\end{matrix} \right\} & (22)\end{matrix}$

Among the real modal coordinates η1l and η*1l introduced in Equation(17) and Equation (18), there are modal coordinates that have a largeeffect on an accuracy of a solution of an original large-scale system(referred to as a dominant mode), and modal coordinates that have asmall effect (referred to as a secondary mode). Subscripts “a” and “b”indicate variables or physical quantities associated with the dominantmodes and the secondary modes, respectively. For example, η1a representsthe dominant modal coordinates and η1b represents the secondary modalcoordinates.

The left and right real modal matrices ^(R)Φ{circumflex over ( )}_(1l)and ^(R)Ψ{circumflex over ( )}_(1l) for the lower-order modes obtainedby Equation (14) are separated into modal matrices ^(R)Φ{circumflex over( )}_(1a) and ^(R)Ψ{circumflex over ( )}_(1a) consisting of n_(a) pairs(2n_(a) pieces) of the dominant modes and modal matrices^(R)Φ{circumflex over ( )}_(1b) and ^(R)Ψ{circumflex over ( )}_(1b)consisting of n_(b) pairs (2n_(b) pieces) of the secondary modes(n_(l)=n_(a)+n_(b)), and they are rearranged as follows. A subscript “#”represents “a” or “b” in the following description.

  R Φ ~ 1 ⁢ l = [   R Φ ~ 1 ⁢ a   R Φ ~ 1 ⁢ b ] ,   R Ψ ~ 1 ⁢ l = [   R Ψ ~1 ⁢ a   R Ψ ~ 1 ⁢ b ] : ( 2 ⁢ n × 2 ⁢ ( n a + n b )   R Φ ~ ? = [ Φ ~ 11 ?  l Φ ~ 1 ?   - l Φ ~ 1 ? Ω 2 ? Φ ~ 22 ? ] ,   R Φ ~ ? = [ Ψ ~ 11 ? l Ψ ~1 ? - l Φ ~ 1 ? Ω 1 2 ? Ψ ~ 22 ? ] : ( 2 ⁢ n × 2 ⁢ n ? ) Φ ~ 11 ? =   R Φ~ ?   l Λ 1 ? -   l Φ ~ ?   R Λ ? , Φ ~ 22 ⁢   ? =   R Φ ~ ?   l Λ ? +  l Φ ~ ?   R Λ ? Ψ ~ 11 ? =   R Ψ ~ 1 ?   l Λ ? -   l Ψ ~ ?   R Λ ? , Ψ ~22 ? =   R Ψ ~ 1 ?   l Λ ? +   l Ψ ~ ?   R Λ ?     R Φ ~ ? =   R Φ ~ ?  l Λ ? ,   l Φ ~ ? =   l Φ ~ ?   l Λ ? ,   R Ψ ~ 1 ? =   R Ψ ~ 1 ?   l Λ? ,   l Ψ ~ 1 ?   =   l Ψ ~ 1 ?     l Λ ? } : ( n × n ? ) } ( 23 )?indicates text missing or illegible when filed

Corresponding to the separation of the real modal matrices, modalcoordinates of Equation (19), Equation (20) and Equation (21) are alsoseparated into the dominant modes and the secondary modes as follows.

1 ξ ¨ 1 ⁢ # + 2 ⁢ Z 1 ⁢ # ⁢ Ω 1 ⁢ # ⁢ ξ . 1 ⁢ # + Ω 1 ⁢ # 2 ⁢ ξ 1 ⁢ # + R 12 ⁢ # ⁢ y¨ 2 + Q 12 ⁢ # ⁢ y . 2 + P 12 ⁢ # ⁢ y 2 =   l Ψ ~ 1 ⁢ # T ⁢ f 1 + Ψ ~ 11 ⁢ # T ⁢f 1 P 12 ⁢ # = - Ω 1 ⁢ # 2 l ⁢ Ψ ~ 1 ⁢ # T ⁢ C 12 + Ψ ~ 11 ⁢ # T ⁢ K 12 , Q 12 ⁢# = - Ω 1 ⁢ # 2 ⁢   l Ψ ~ 1 ⁢ # 2 ⁢ M 12 + Ψ ~ 1 ⁢ # T ⁢ M 12 + Ψ ~ 22 ⁢ # T ⁢ C12 +   l Ψ ~ 1 ⁢ # T ⁢ K 12 R 12 ⁢ # = Ψ ~ 22 ⁢ # T ⁢ M 12 , P 12 ⁢ # , Q 12 ⁢# , R 12 ⁢ # : ( n ? × m ) } ( 24 ) $\begin{matrix}\left. \begin{matrix}{{{\overset{¨}{\xi}}_{1\#}^{*} + {2Z_{1\#}\Omega_{1\#}{\overset{.}{\xi}}_{1\#}^{*}} + {\Omega_{1\#}^{2}\xi_{1\#}^{*}} + {R_{12\#}^{*}{\overset{¨}{y}}_{2}} + {Q_{12\#}^{*}{\overset{.}{y}}_{2}} + {P_{12\#}^{*}y_{2}}} = 0} \\\begin{matrix}{{P_{12\#}^{*} = {{{- \Omega_{1\#}^{2}}{\,^{l}{\overset{\sim}{\Omega}}_{1\#}^{T}}C_{21}^{T}} + {{\overset{\sim}{\Phi}}_{11\#}^{T}K_{21}^{T}}}},} \\{Q_{12\#}^{*} = {{{- \Omega_{1\#}^{2}}{\,^{l}{\overset{\sim}{\Phi}}_{1\#}^{T}}M_{21}^{T}} + {{\overset{\sim}{\Phi}}_{22\#}^{T}C_{21}^{T}} + {{\,^{l}{\overset{\sim}{\Phi}}_{1\#}^{T}}K_{21}^{T}}}}\end{matrix} \\{{R_{12\#}^{*} = {{\overset{\sim}{\Phi}}_{22\#}^{T}M_{21}^{T}}},P_{12\#}^{*},Q_{12\#}^{*},{R_{12\#}^{*}:\left( {n_{\#} \times m} \right)}}\end{matrix} \right\} & (25)\end{matrix}$ $\begin{matrix}\left. \begin{matrix}{{\left( {M_{22} + {\overset{\sim}{M}}_{22h}} \right){\overset{¨}{y}}_{2}} + {\left( {C_{22} + {\overset{\sim}{C}}_{22h}} \right){\overset{.}{y}}_{2}} + {\left( {K_{22} + {\overset{\sim}{K}}_{22h}} \right)y_{2}} +} \\{{R_{12a}^{*}\text{?}{\overset{¨}{\xi}}_{1a}} + {Q_{12a}^{*T}{\overset{.}{\xi}}_{1a}} + {P_{12a}^{T}\xi_{1a}} + {R_{12b}^{*T}{\overset{¨}{\xi}}_{1b}} + {Q_{12b}^{*T}{\overset{.}{\xi}}_{1b}} + {P_{12b}^{*T}\xi_{1b}} +} \\{n_{2} = {f_{2} + {q_{21h}{\overset{.}{f}}_{1}} + {p_{21h}f_{1}}}}\end{matrix} \right\} & (26)\end{matrix}$ ?indicates text missing or illegible when filed

Since an effect of nonlinearity is small for the secondary modes, then,ignoring this effect and assuming that y₂, ξ_(1b) and ξ*_(1b) vibrateharmonically with the fundamental period 2π/ω yields y{umlaut over( )}2=−ω²y₂, ξ{umlaut over ( )}_(1b)=−ω²ξ_(1b), ξ*{umlaut over( )}_(1b)=−ω²ξ*_(1b). Furthermore, considering that f₁=−ω²f₁ holds forthe harmonic external force, from Equations (24) and (25), anapproximate solution of the secondary modes is obtained as follows.

$\begin{matrix}\left. \begin{matrix}{\begin{bmatrix}\xi_{1b} \\{\overset{.}{\xi}}_{1b}\end{bmatrix} = {{\begin{bmatrix}U_{12b} & V_{12b} \\{{- \omega^{2}}V_{12b}} & U_{12b}\end{bmatrix}\begin{bmatrix}y_{2} \\{\overset{.}{y}}_{2}\end{bmatrix}} + {\begin{bmatrix}U_{11b} & V_{11b} \\{{- \omega^{2}}V_{11b}} & U_{11b}\end{bmatrix}\begin{bmatrix}f_{1} \\{\overset{.}{f}}_{1}\end{bmatrix}}}} \\{U_{12b} = {{- \Delta_{1b}^{- 1}}\left\{ {{\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right)\left( {P_{12b} - {\omega^{2}R_{12b}}} \right)} - {2\omega^{2}{\,^{R}\Lambda_{1b}}Q_{12b}}} \right\}}} \\{V_{12b} = {{- \Delta_{1b}^{- 1}}\left\{ {{2{\,^{R}\Lambda_{1b}}\left( {P_{12b} - {\omega^{2}R_{12b}}} \right)} + {\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right)Q_{12b}}} \right\}}} \\{U_{11b} = {\Delta_{1b}^{- 1}\left\{ {{\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right){\overset{\sim}{\Psi}}_{11b}^{T}} - {2\omega^{2}{\,^{R}\Lambda_{1b}}{\,^{l}{\overset{\sim}{\Psi}}_{1b}^{T}}}} \right\}}} \\{V_{11b} = {\Delta_{1b}^{- 1}\left\{ {{2{\,^{R}\Lambda_{1b}}{\overset{\sim}{\Psi}}_{11b}^{T}} + {\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right){\,^{l}{\overset{\sim}{\Psi}}_{1b}^{T}}}} \right.}} \\{\Delta_{1b} = {\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right)^{2} + \left( {2\omega{\,^{R}\Lambda_{1b}}} \right)^{2}}}\end{matrix} \right\} & (27)\end{matrix}$ $\begin{matrix}\left. \begin{matrix}{\begin{bmatrix}\xi_{1b}^{*} \\{\overset{\sim}{\xi}}_{1b}^{*}\end{bmatrix} = {\begin{bmatrix}U_{12b}^{*} & V_{12b}^{*} \\{{- \omega^{2}}V_{12b}^{*}} & U_{12b}^{*}\end{bmatrix}\begin{bmatrix}y_{2} \\{\overset{.}{y}}_{2}\end{bmatrix}}} \\{U_{12b}^{*} = {{- \Delta_{1b}^{- 1}}\left\{ {{\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right)\left( {P_{12b}^{*} - {\omega^{2}R_{12b}^{*}}} \right)} - {2\omega^{2}{\,^{R}\Lambda_{1b}}Q_{12b}^{*}}} \right\}}} \\{V_{12b}^{*} = {{- \Delta_{1b}^{- 1}}\left\{ {{2{\,^{R}\Lambda_{1b}}\left( {P_{12b}^{*} - {\omega^{2}R_{12b}^{*}}} \right)} + {\left( {\Omega_{1b}^{2} - {\omega^{2}I_{n_{b}}}} \right)Q_{12b}^{*}}} \right\}}}\end{matrix} \right\} & (28)\end{matrix}$

By substituting the approximate solution of Equation (27) into thesecondary modal coordinates ξ_(1b), ξ{dot over ( )}_(1b), ξ{umlaut over( )}_(1b) in Equation (26) and correcting and eliminating the secondarymodes, the following equation is obtained (the result of Equation (28)is used in this derivation process).

$\begin{matrix}\left. \begin{matrix}{{\left( {M_{22} + {\overset{\sim}{M}}_{22b} + {\overset{\sim}{M}}_{22b}} \right){\overset{¨}{y}}_{2}} + {\left( {C_{22} + {\overset{\sim}{C}}_{22b} + {\overset{\sim}{C}}_{22b}} \right){\overset{.}{y}}_{2}} +} \\{{\left( {K_{22} + {\overset{\sim}{K}}_{22b} + {\overset{\sim}{K}}_{22b}} \right)y_{2}} + {R_{12a}^{*T}{\overset{¨}{\xi}}_{1a}} + {Q_{12a}^{*T}{\overset{.}{\xi}}_{1a}} +} \\{{{P_{12a}^{*T}\xi_{1a}} + n_{2}} = {f_{2} + {\left( {q_{21b} + q_{21b}} \right){\overset{.}{f}}_{1}} + {\left( {p_{21b} + p_{21b}} \right)f_{1}}}}\end{matrix} \right\} & (29)\end{matrix}$

Where, M{tilde over ( )}_(2b), C{tilde over ( )}_(2b), K{tilde over( )}₂₂, p_(21b) and q_(2lb) represent correction terms of the secondarymodes and are given by the following equations.

$\begin{matrix}\left. \begin{matrix}{{\overset{\sim}{M}}_{22b} = {{M_{21}{\,^{l}{\overset{\sim}{\Phi}}_{1l}}\Omega_{1l}^{2}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}M_{12}} + {U_{12b}^{*T}U_{12b}} - {V_{12b}^{*T}\Omega_{1b}^{2}V_{12b}} +}} \\{{U_{12b}^{*T}R_{12b}} + {R_{12b}^{*T}U_{12b}}} \\{{\overset{\sim}{C}}_{22b} = {{{- P_{12l}^{*T}}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}M_{12}} + {M_{21}^{l}{\overset{\sim}{\Phi}}_{1l}\Omega_{1l}^{2l}{\overset{\sim}{\Psi}}_{1l}^{T}C_{12}} + {R_{12l}^{*T}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}K_{12}} +}} \\{{\left( {P_{12b}^{*T} - {\omega^{2}R_{12b}^{*T}}} \right)V_{12b}} + {Q_{12b}^{*T}U_{12b}}} \\{{\overset{\sim}{K}}_{22b} = {{{- P_{12l}^{*T}}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}C_{12}} + {\left( {Q_{12l}^{*T} + {{M_{21}}^{l}{\overset{\sim}{\Phi}}_{1l}\Omega_{1l}^{2}}} \right){\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}K_{12}} +}} \\{{\left( {P_{12b}^{*T} - {\omega^{2}R_{12b}^{*T}}} \right)U_{12b}} - {\omega^{2}Q_{12b}^{*T}V_{12b}} + {\omega^{2}\left( {{U_{12b}^{*T}U_{12b}} -} \right.}} \\\left. {{V_{12b}^{*T}\Omega_{1b}^{2}V_{12b}} + {U_{12b}^{*T}R_{12b}} + {R_{12b}^{*T}U_{12b}}} \right) \\{p_{21b} = {{\left( {Q_{12l}^{*T} + {M_{21}{\,^{l}{\overset{\sim}{\Phi}}_{1l}}\Omega_{1l}^{2}}} \right){\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}} - {\left( {P_{12b}^{*T} - {\omega^{2}R_{12b}^{*T}}} \right)U_{11b}} + {\omega^{2}Q_{12b}^{*T}V_{11b}}}} \\{q_{21b} = {{R_{12l}^{*T}{\,^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}} - {\left( {P_{12b}^{*T} - {\omega^{2}R_{12b}^{*T}}} \right)V_{11b}} - {Q_{12b}^{*T}U_{11b}}}}\end{matrix} \right\} & (30)\end{matrix}$

Furthermore, from the equation with the subscript W “#” in Equation (24)replaced with “a” and Equation (29), an equation whose dimensionality isreduced to (n_(a)+m) dimension is obtained as follows.

$\begin{matrix}\left. \begin{matrix}{{{M\overset{¨}{\xi}} + {C\overset{.}{\xi}} + {K\xi} + n} = f} \\{{\xi = \begin{bmatrix}\xi_{1a} \\y_{2}\end{bmatrix}},{M = \begin{bmatrix}I_{n_{i}} & R_{12a} \\R_{12a}^{*T} & {M_{22} + {\overset{\sim}{M}}_{22}}\end{bmatrix}},} \\{C = \begin{bmatrix}{2Z_{1a}\Omega_{1a}} & Q_{12a} \\Q_{12a}^{*T} & {C_{22} + {\overset{\sim}{C}}_{22}}\end{bmatrix}} \\{{K = \begin{bmatrix}\Omega_{1a}^{2} & P_{12a} \\P_{12a}^{*T} & {K_{22} + {\overset{\sim}{K}}_{22}}\end{bmatrix}},{n = \begin{bmatrix}0 \\n_{2}\end{bmatrix}},{f = \begin{bmatrix}{{{\,^{l}{\overset{\sim}{\Psi}}_{1a}^{T}}{\overset{.}{f}}_{1}} + {{\overset{\sim}{\Psi}}_{11a}^{T}f_{1}}} \\{f_{2} + {q_{21}{\overset{.}{f}}_{1}} + {p_{21}f_{1}}}\end{bmatrix}}} \\{{{\overset{\sim}{M}}_{22} = {{\overset{\sim}{M}}_{22b} + {\overset{\sim}{M}}_{22h}}},{{\overset{\sim}{C}}_{22} = {{\overset{\sim}{C}}_{22b} + {\overset{\sim}{C}}_{22h}}},{{\overset{\sim}{K}}_{22} = {{\overset{\sim}{K}}_{22b} + {\overset{\sim}{K}}_{22b}}}} \\{{q_{21} = {q_{21b} + q_{21h}}},{p_{21} = {p_{21b} + p_{21h}}}}\end{matrix} \right\} & (31)\end{matrix}$

The model represented by Equation (31) is called the dimension reducedmodel. The approximate solution is calculated using Equation (31).

When approximate solutions (ξ_(1a), ξ{dot over ( )}_(1a), y₂, y{dot over( )}₂ are obtained from Equation (31) of the dimension reduced model,the approximate solutions for linear state variables are converted intophysical coordinates when necessary. The relational equations requiredfor this procedure are as follows.

$\begin{matrix}{\begin{bmatrix}y_{1} \\{\overset{.}{y}}_{1}\end{bmatrix} = {\begin{bmatrix}y_{1a} \\{\overset{.}{y}}_{1a}\end{bmatrix} + \begin{bmatrix}y_{1b} \\{\overset{.}{y}}_{1b}\end{bmatrix} + \begin{bmatrix}y_{1h} \\{\overset{.}{y}}_{1h}\end{bmatrix}}} & (32)\end{matrix}$ Where, $\begin{matrix}{\begin{bmatrix}y_{1a} \\{\overset{.}{y}}_{1a}\end{bmatrix} = {\begin{bmatrix}{\overset{\sim}{\Phi}}_{11a} & {\,^{l}{\overset{\sim}{\Phi}}_{1a}} \\{{\,^{l}{\overset{\sim}{\Phi}}_{1a}}\Omega_{1a}^{2}} & {\overset{\sim}{\Phi}}_{22a}\end{bmatrix}\left\{ {\begin{bmatrix}\xi_{1a} \\{\overset{.}{\xi}}_{1a}\end{bmatrix} + {\begin{bmatrix}{{- {\,^{l}{\overset{\sim}{\Psi}}_{1a}^{T}}}C_{12}} & {{- {\,^{l}{\overset{\sim}{\Psi}}_{1a}^{T}}}M_{12}} \\{{\,^{l}{\overset{\sim}{\Psi}}_{1a}^{T}}K_{12}} & 0\end{bmatrix}\begin{bmatrix}y_{2} \\{\overset{\sim}{y}}_{2}\end{bmatrix}} + {\begin{bmatrix}0 & 0 \\{- {\,^{l}{\overset{\sim}{\Psi}}_{1a}^{T}}} & 0\end{bmatrix}\begin{bmatrix}{\overset{.}{f}}_{1} \\f_{1}\end{bmatrix}}} \right\}}} & (33)\end{matrix}$ $\begin{matrix}\left. \begin{matrix}{\begin{bmatrix}y_{1b} \\{\overset{.}{y}}_{1b}\end{bmatrix} = {\begin{bmatrix}{\overset{\sim}{\Phi}}_{11b} & {\,^{l}{\overset{\sim}{\Phi}}_{1b}} \\{{- {\,^{l}{\overset{\sim}{\Phi}}_{1b}}}\Omega_{1b}^{2}} & {\overset{\sim}{\Phi}}_{22b}\end{bmatrix}\left\{ \begin{bmatrix}{U_{12b} - {{\,^{l}{\overset{\sim}{\Psi}}_{1b}^{T}}C_{12}}} & {V_{12b} - {{\,^{l}{\overset{\sim}{\Psi}}_{1b}^{T}}M_{12}}} \\{{{- \omega^{2}}V_{12b}} + {{\,^{l}{\overset{\sim}{\Psi}}_{1b}^{T}}K_{12}}} & U_{12b}\end{bmatrix} \right.}} \\\left. {{{\begin{bmatrix}y_{2} \\{\overset{.}{y}}_{2}\end{bmatrix}++}\begin{bmatrix}U_{11b} & V_{11b} \\{{{- \omega^{2}}V_{11b}} - {\,^{l}{\overset{\sim}{\Psi}}_{1b}^{T}}} & U_{11b}\end{bmatrix}}\begin{bmatrix}f_{1} \\{\overset{.}{f}}_{1}\end{bmatrix}} \right\}\end{matrix} \right\} & (34)\end{matrix}$ $\begin{matrix}{\begin{bmatrix}y_{1h} \\{\overset{.}{y}}_{1h}\end{bmatrix} = {{- {\begin{bmatrix}{\Gamma_{2h}K_{12}} & {{\Gamma_{2h}C_{12}} + {\Gamma_{3h}K_{12}}} \\{{- \Gamma_{1h}}K_{12}} & {{{- \Gamma_{1h}}C_{12}} + {\Gamma_{2h}K_{12}}}\end{bmatrix}\begin{bmatrix}y_{2} \\{\overset{.}{y}}_{2}\end{bmatrix}}} + {\begin{bmatrix}\Gamma_{2h} & \Gamma_{3h} \\{- \Gamma_{1b}} & \Gamma_{2h}\end{bmatrix}\begin{bmatrix}f_{1} \\{\overset{.}{f}}_{1}\end{bmatrix}}}} & (35)\end{matrix}$

y_(1a), y{dot over ( )}_(1b) represent physical coordinates based on thesecondary modes, and y_(1h), y{dot over ( )}_(1h) represent physicalcoordinates based on the higher-order modes.

When full content of a complicated and diverse frequency response isclarified while gradually changing an angular frequency ω of an externalforce in a nonlinear system, an approximate solution obtained for acertain ω is often used as an initial value for ω±Δω in iterativeapproximation calculation. Therefore, ξ_(1a) to be used for ω±Δω can beselected by using information included in the approximate solutionobtained for a certain ω.

The relational equations required for this procedure are as follows.

Now, suppose that an approximate solution of the dimension reduced modelconsisting of ξ_(1a), ξ{dot over ( )}_(1a), y₂ and y{dot over ( )}₂ isobtained for a certain ω. From this approximate solution y₂, y{dot over( )}₂, an approximate solution of lower-order modal coordinates ξ_(1l)can be calculated by the following equation.

$\begin{matrix}\left. \begin{matrix}{\begin{bmatrix}\xi_{1l} \\{\overset{.}{\xi}}_{1l}\end{bmatrix} = {{\begin{bmatrix}U_{12l} & V_{12l} \\{{- \omega^{2}}V_{12l}} & U_{12l}\end{bmatrix}\begin{bmatrix}y_{2} \\{\overset{.}{y}}_{2}\end{bmatrix}} + {\begin{bmatrix}U_{11l} & V_{11l} \\{{- \omega^{2}}V_{11l}} & U_{11l}\end{bmatrix}\begin{bmatrix}f_{1} \\{\overset{.}{f}}_{1}\end{bmatrix}}}} \\{U_{12l} = {{- \Delta_{1l}^{- 1}}\left\{ {{\left( {\Omega_{1l}^{2} - {\omega^{2}I_{n_{l}}}} \right)\left( {P_{12l} - {\omega^{2}R_{12l}}} \right)} - {2\omega^{2R}\Lambda_{1l}Q_{12l}}} \right\}}} \\{V_{12l} = {{- \Delta_{1l}^{- 1}}\left\{ {{2^{R}{\Lambda_{1l}\left( {P_{12l} - {\omega^{2}R_{12l}}} \right)}} + {\left( {\Omega_{1l}^{2} - {\omega^{2}I_{n_{l}}}} \right)Q_{12l}}} \right\}}} \\{U_{11l} = {\Delta_{1l}^{- 1}\left\{ \left( {\Omega_{1l}^{2}\left\{ {{\left( {\Omega_{1l}^{2} - {\omega^{2}I_{n_{l}}}} \right){\overset{\sim}{\Psi}}_{11l}^{T}} - {2\omega^{2R}\Lambda_{1l}^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}} \right\}} \right. \right.}} \\{V_{11l} = {\Delta_{1l}^{- 1}\left\{ {{2^{R}\Lambda_{1l}{\overset{\sim}{\Psi}}_{11l}^{T}} + {\left( {\Omega_{1l}^{2} - {\omega^{2}I_{n_{l}}}} \right)^{l}{\overset{\sim}{\Psi}}_{1l}^{T}}} \right\}}} \\{\Delta_{1l} = {\left( {\Omega_{1l}^{2} - {\omega^{2}I_{n_{l}}}} \right)^{2} + \left( {2\omega^{R}\Lambda_{1l}} \right)^{2}}}\end{matrix} \right\} & (36)\end{matrix}$

An amplitude for each mode of this approximate solution can be definedby the following equation.

$\begin{matrix}{{\Delta\xi_{1,p}} = {{\xi_{1,p}} = {\sqrt{\frac{\omega}{\pi}{\int}_{0}^{2\pi/\omega}\left\{ {\xi_{1,p}(t)} \right\}^{2}{dt}}\left( {{p = 1},\ldots,n_{l}} \right)}}} & (37)\end{matrix}$

In the calculation of ω±Δω, only a mode whose modal amplitude ∥ξ_(1,p)∥is greater than a set threshold value δ is selected as ξ_(1a). If it isassumed that the number of selected dominant modes is n_(a) and thenumber of remaining secondary modes is n_(b), n_(l)=n_(a)+n_(b). Anupper limit value n_(a,max) of n_(a) may be defined depending on a case.

With this method, it is possible to efficiently calculate a frequencyresponse while appropriately selecting dominant modes as few as possibleaccording to characteristics of a solution. In this method, the problemis how to set the threshold value δ for the modal amplitude and theupper limit value n_(a,max) for the number of dominant modes. Inpractical use, by obtaining multiple frequency responses while graduallydecreasing δ and gradually increasing n_(a,max), and when the analyticalresults no longer change, a highly accurate solution is considered tohave been obtained. Since the dimension reduced model can be calculatedfairly fast, it is considered that a cost for a plurality ofcalculations is not a significant issue in practice. Note that a casewhen the analytical results do not change includes not only a case wherethe analytical results do not change at all, but also a case where adifference between the analytical results is less than a predeterminedreference.

As for the method of setting the number of dominant modes n_(a) at afirst ω at which the calculation is started, for example, n_(a)=n₁(thenumber of lower-order modes) may be set.

A procedure for calculating the frequency response is briefly summarizedbelow (refer to FIG. 1 and FIG. 2 ).

-   -   (Step 1) Preprocessing    -   (Step 1-1) Calculation of the complex eigenvalues and the        complex eigenvectors of the lower-order modes for the linear        state variables [Equation (8), Equation (10), Equation (12),        Equation (13)] and realization of the complex eigenvectors        [Equation (14)]    -   (Step 1-2) Introduction of the real modal coordinates for the        physical coordinates of the lower-order modes [Equation (17),        Equation (18)] and derivation of the real modal equation        [Equation (19), Equation (20)]    -   (Step 1-3) Correction of the effect of the higher-order modes        using the lower-order modes for the equation for the nonlinear        state variables [Equation (21)]    -   (Step 2) Selection of the dominant mode [Equation (36), Equation        (37)]    -   (Step 3) Derivation of the dimension reduced model    -   (Step 3-1) Separation of the dominant mode and the secondary        mode [Equation (23), Equation (24), Equation (25), Equation        (26)]    -   (Step 3-2) Calculation of the approximate solution of the        secondary mode [Equation (27), Equation (28)]    -   (Step 3-3) Derivation of the dimension reduced model [Equation        (29), Equation (30), Equation (31)]    -   (Step 4) Calculation of the periodic solution of the dimension        reduced model [Equation (31)].    -   (Step 5) Determination of calculation end or continuation.

The calculation will end when the angular frequency ω reaches an upperor lower limit of an analysis range. When it does not reach the limit,ω±Δω is replaced with ω. The processing returns to Step 1 when M₁₁, C₁₁,and K₁₁ are the functions of ω, and returns to Step 2 to continue thecalculation when they are constants.

[1. Analytical model] To confirm effectiveness of the dimensionalityreduction method using the new type of complex modal analysis,particularly effectiveness of the elimination method of higher-ordermodes and the selection method of dominant modes, the dimensionalityreduction method is applied to an in-plane bending vibration of astraight beam structure. FIG. 3 shows a schematic diagram of ananalytical model, and Table 1 shows parameters used for the calculation.

TABLE 1 Beam Total length [m]  1.5 Outer diameter [mm] 30 Inner diameter[mm] 25 Distributed Translational direction 50 damping of [kg/(m s)]beam Rotational direction  1.0 [kg/s] Restoring force of support portionBoth ends [N] Position of 0.5 m from left end [N] 10³{dot over (y)} +10⁷y + 10¹²y³ $\begin{matrix}{{\left. \begin{matrix}{0:} & {{❘y❘} \leq \delta} \\{10^{7}\left( {y - \Delta} \right):} & {{❘y❘} > \delta}\end{matrix} \right\}\delta} =} \\{{2.5{mm}},{\Delta = {{\delta \cdot {sign}}(y)}}}\end{matrix}$ External force Position of 1 m 0.001 ω²sinωt from left end[N]

In this model 10, both ends of a beam 1 having a uniform hollow circularsection of 1.5 m in length, 30 mm in outer diameter, and 25 mm in innerdiameter are supported by nonlinear elements 2 (viscous dampers andcubic springs of hardening type), and a piecewise linear spring 3 (gap)is disposed at a position of 0.5 m from the left end. In addition, aharmonic external force F of a centrifugal force type is applied on aposition of 1 m from the left end so that a vibration with a largeamplitude generates and the effect of nonlinearity strongly appear evenin a region where a frequency of the external force is relatively high.

Physical properties of beam 1 are assumed to be steel material. The beam1 has been equally divided into 30 small elements, and the mass andstiffness matrices have been derived by the finite element method. Inthis case, the degrees of freedom are n=59 for the linear statevariables, m=3 for the nonlinear state variables, and N=62 for the totaldegree of freedom

The degrees of freedom of this analytical model are not large enough tobe called a large-scale system, but they are close to a limit at whichthe solution can be obtained without applying the dimensionalityreduction method. The damping matrix is assumed to be a non-proportionaldamping based on distributed damping so that an advantage of the newtype of complex modal analysis appears. Table 2 shows the first-order totenth order undamped natural frequencies of the linearized system withthe coefficient of the cubic term of the nonlinear spring in thisanalytical model set to zero.

TABLE 2 Order number Natural frequency [Hz] 1 34.81 2 138.2 3 306.9 4535.1 5 812.7 6 1125 7 1459 8 1820 9 2235 10 2728

FIGS. 4 to 9 show the frequency response diagrams of this analyticalmodel. FIGS. 4, 6, and 8 show the first-order to fifth-order primaryresonance regions, and FIGS. 5, 7, and 9 show the enlarged views of thefirst-order and second-order primary resonance regions. FIGS. 4 and 5show the results when the dimensionality reduction method is not applied(referred to as the full model), and FIGS. 6 to 9 show results when thedimensionality reduction method is applied (referred to as the dimensionreduced model).

The dimension reduced model is constructed by selecting a fixed number(n_(a)=5, 15) of dominant modes in the order that the natural frequencyof the constraint mode is closer to the frequency of the external force,without both applying the elimination method of higher-order modes andthe selection method of dominant modes. In the results of the dimensionreduced model, the numbers of dominant modes (numerical values aredenoted on the right axis) are indicated by a filled V marks, but theyare constant in FIGS. 6 to 9 , so that the results are presented as anapparently thick solid line.

All response curves in the frequency response denote a maximumhalf-amplitude at the point where the gap is disposed and bend when theamplitudes exceed a gap width (0.0025 m). In the figures, a solid linerepresents a stable solution, a dashed line represents an unstablesolution, an ◯ mark represents a saddle-node bifurcation point, a □ markrepresents a pitchfork bifurcation point, and a Δ mark represents a hopfbifurcation point.

A result of comparison between FIGS. 6 and 7 when n_(a)=5 and FIGS. 4and 5 of the full model shows a generally good consistency therebetween,but the shapes of the response curves, the peak amplitude values, andthe results of stability analysis are different in a region where theeffect of a gap that is strong nonlinearity appears. As shown in FIGS. 8and 9 , when the number of dominant modes is increased to n_(a)=15, theresults consistent with those of the full model are obtained.

[2. Application of elimination method of higher-order modes] Theelimination method of higher-order modes is applied to this analyticalmodel, and the effect on an analytical result of the number oflower-order modes n_(l) used for the analysis is investigated. FIGS. 10and 11 show results when the number of dominant modes is fixed atn_(a)=5 and the number of lower-order modes is n_(l)=5, 10. Compared toFIG. 6 (n_(l)=n=59) when the elimination method of higher-order modes isnot applied, almost the same response is obtained in the first-order andsecond-order primary resonance regions with n_(l)=5 in FIG. 10 .However, the difference increases as the frequency increases. Ifn_(l)=10 in FIG. 11 , it can be seen that the almost consistent resultsare obtained in the analytical frequency range.

In this manner, the required number of lower-order modes changesdepending on the range of frequency for performing the analysis, and thesolution can be obtained with high accuracy up to a higher frequencyrange by increasing the number of lower-order modes. In this analyticalmodel, since the effect of higher-order mode can be estimated with highaccuracy with only 10 out of 59, that is a small number of lower-ordermodes, this property is expected to be extremely effective for ananalysis of large-scale systems.

To examine the effect of the number of lower-order modes on the naturalfrequencies, the first-order to fifth-order natural frequencies in theanalytical frequency range are shown in FIGS. 12 to 16 . These naturalfrequencies were obtained from the dimension reduced model (n_(l)=n_(a))derived by applying only the elimination method of higher-order modes.The dots represent the natural frequencies of the dimension reducedmodel, and the dashed line represents the natural frequency of the fullmodel.

From the comparison of FIGS. 12 to 16 , the lower order naturalfrequencies have better accuracy and the higher the order, the lower theaccuracy. It can also be seen that the accuracy improves as the numberof lower-order modes increases. Since the accuracy of the naturalfrequencies and the accuracy of the frequency response analyticalresults show similar trends with respect to the number of lower-ordermodes, it is considered that the accuracy of the natural frequencies canbe used as a guide for setting the number of lower-order modes n_(l)from the viewpoint of calculation cost.

[3. Application of selection method of dominant modes] The selectionmethod of dominant modes is applied to the analytical model. Theanalysis range is limited to the second-order primary resonance region.The elimination method of higher-order modes is also applied to thedimensionality reduction method, and the number of lower-order modes isset as n_(l)=20. FIG. 17 is a result of the full model, and FIGS. 18 to20 are results of applying the selection method of dominant modes withthe respective threshold values set as δ=5×10⁻⁴, 1×10⁻⁴, and 1×10⁻⁵,respectively. It can be seen that the number of dominant modes varieswith the frequency of the external force. It can be seen that when thethreshold value is gradually decreased from FIG. 18 to FIG. 20 , thenumber of modes increases and the modes required for the analysis areappropriately selected, and thereby the solution of the full model isapproached. At δ=1×10⁻⁵ in FIG. 20 , the results consistent with thoseof the full model are obtained.

It is confirmed that, at any threshold value, the number of dominantmodes is relatively small in a region where the amplitude of thesolution is small and the effect of nonlinearity does not appear, thenumber of dominant modes increases as an amplitude of a solutionincreases, and the analysis can be performed efficiently by changing thedominant modes according to the characteristics of the solution.

[4. Model with large degree of freedom] The present method is applied tothe analytical model with 300 divisions to verify the effectiveness. Inthis case, the degrees of freedom are n=599 for the linear statevariables, m=3 for the nonlinear state variables, and N=602 for thetotal degree of freedom. With this degree of freedom, the analysis forthe full model becomes difficult. The analysis is performed separatelyin the first-order to second-order primary resonance regions and thefourth-order to fifth-order primary resonance regions.

[4.1. First-order and second-order primary resonance regions] Theanalytical range of the frequency of the external force is set to 20 Hzto 200 Hz based on the natural frequencies and the results of Chapter 1[1. analytical model]. The number of lower-order modes when theelimination method of higher-order modes is applied is set to n_(l)=15based on a sixth-order natural frequency, which is about five times themaximum frequency of 200 Hz.

FIGS. 21 to 23 show the results of applying the selection method ofdominant modes while respective threshold values are gradually decreasedto δ=7×10⁻⁴, 1×10⁻⁴, and 2×10⁻⁵. A comparison between FIG. 21 and FIG.22 shows that there is a difference in the responses of the first-orderand second-order resonances. From a comparison between FIG. 22 and FIG.23 , it can be seen that a solution near the second-order peak isslightly different although the change in response is small

To confirm that the result of FIG. 23 is a highly accurate result, ananalysis of the dimension reduced model is performed while n_(a) isgradually increased as in Chapter 1 without applying both theelimination method of higher-order modes and the selection method ofdominant modes, and a result when there was no change in the responsewas regarded as a correct solution for the analytical model. The resultis shown in FIG. 24 , where the number of lower-order modes is n_(l)=599and the number of dominant modes is n_(a)=20. It can be seen that thisresult shows very good consistency with FIG. 23 .

[4.2. Fourth-order and fifth-order primary resonance regions] Theanalytical frequency range is set from 450 Hz to 1000 Hz, and the numberof lower-order modes is set to n_(l)=30 based on the tenth-order naturalfrequency, which is about three times the maximum frequency of 1000 Hz.As in the previous section, the selection method of dominant modes isapplied while the threshold value is gradually decreased. FIGS. 25 to 27show the results when the threshold values are set to δ=1×10⁻³, 5×10⁻⁴,and 2×10⁻⁴, respectively.

From the comparison of FIGS. 25 to 27 in the same manner as in theprevious section, it can be seen that the change in response of thedimension reduced model is getting smaller. In addition, from thecomparison between FIG. 27 and the results of n_(l)=599 and n_(a)=20 inFIG. 28 , even in this frequency range, the results that are almost thesame as those of the full model can be efficiently calculated by thedimensionality reduction method. However, the number of lower-ordermodes and the optimum threshold value differ depending on the frequencyrange for performing the analysis.

As described above, the effectiveness of the present invention wasconfirmed through specific numerical calculations for this analyticalmodel. Note that the analytical method described above may beimplemented by an analysis device 100 as shown in FIG. 29 .

The analysis device 100 is realized by a device such as a personalcomputer, a server, or an industrial computer. The analysis device 100includes a processing unit 110, a storage unit 120, an input unit 130,and an output unit 140. The processing unit 110 applies the new type ofcomplex modal analysis to the equation for linear state variables toconvert it into the real modal equations for lower-order modes andcorrects and eliminates the effect of higher-order modes of the linearstate variables from the equation for nonlinear state variable. Theprocessing unit 110 selects the dominant modes, which have the largeeffect on the solution of the original large-scale system, from the realmodal equation, and, in relation to the secondary mode, which have thesmall effect, eliminates the modes thereof by incorporating the effectto the equation for nonlinear state variables as correction termsobtained from an approximate solution of the real modal equation forlower-order modes, and derives the dimension reduced model. Theprocessing unit 110 calculates the frequency response by using thedimension reduced model. The processing unit 110 is realized by, forexample, a hardware processor such as a central processing unit (CPU)executing a program (software) stored in the storage unit 120. Inaddition, some or all of these functional units may be realized byhardware (a circuit unit; including circuitry) such as large scaleintegration (LSI), an application specific integrated circuit (ASIC), anfield-programmable gate array (FPGA), and a graphics processing unit(GPU), or may be realized by software and hardware in cooperation. Theprogram may be stored in advance in a storage device (a storage deviceincluding a non-transitory storage medium) such as a hard disk drive(HDD) or flash memory or may be stored in a removable storage medium(non-transitory storage medium) such as a DVD or CD-ROM and installed bythe storage medium being attached to a drive device. The input unit 130is realized by, for example, a keyboard, a mouse, a touch panel, or thelike. The output unit 140 is realized by, for example, a display, aprinter, a touch panel, or the like. Men an analysis method is realized,information that needs to be set, such as a threshold value, may bestored in advance in the storage unit 120 or may be input by aresearcher from the input unit 130. The analytical result may be outputto the output unit 140.

1. A vibration analysis method for analyzing vibrations of a large-scalesystem with local strong nonlinearities comprising: a process (1) ofapplying the new type of complex modal analysis to an equation forlinear state variables to convert the equation to a real modal equationfor lower-order modes, and correcting an effect of higher-order modes ofthe linear state variables and eliminating the modes from an equationfor nonlinear state variables; a process (2) of selecting dominantmodes, which have a large effect on a solution of an originallarge-scale system, from the real modal equation for lower-order modes,and, in relation to secondary modes, which have a small effect,eliminating the modes by incorporating the effect to the equation fornonlinear state variables as a correction term obtained from anapproximate solution of the real modal equation for lower-order modes,and deriving the dimension reduced model; and a process (3) ofcalculating a frequency response by using the dimension reduced model.2. The vibration analysis method according to claim 1, wherein, inrelation to an angular frequency ω, an angular frequency ω±Δω isreplaced with ω and the process (1) to the process (3) are repeateduntil the angular frequency ω reaches an upper or lower limit of ananalysis range.
 3. The vibration analysis method according to claim 2,wherein, when an angular frequency ω±Δω is replaced with ω and theprocess (1) to the process (3) are repeated, the dominant modes and thesecondary modes at an angular frequency ω±Δω are separated on the basisof the approximate solution obtained for the dimension reduced model atan angular frequency ω.
 4. The vibration analysis method according toclaim 3, wherein the dominant modes and the secondary modes at anangular frequency ω±Δω are separated on the basis of a relationshipbetween a modal amplitude obtained from the dimension reduced model atan angular frequency ω and a predetermined threshold value.
 5. Thevibration analysis method according to claim 1, wherein, in the new typeof complex modal analysis, calculation of complex eigenvalues andcomplex eigenvectors for lower-order modes of linear state variables andrealization of the complex eigenvectors are performed, and introductionof real modal coordinates for physical coordinates of lower-order modesand derivation of real modal equations are performed.
 6. A program whichcauses a computer as a device for analyzing vibrations of a large-scalesystem with local strong nonlinearities to function as: a firstprocessing unit configured to apply the new type of complex modalanalysis to an equation for a linear state variable to convert theequation to a real modal equation for lower-order modes and correct aneffect of higher-order modes of the linear state variables and eliminatethe modes from an equation for a non-linear state variable, a secondprocessing unit configured to select dominant modes, which have a largeeffect on a solution of an original large-scale system, from the realmodal equation for lower-order modes, and, in relation to secondarymodes, which have a small effect, eliminate the modes by incorporatingthe effect to the equation for nonlinear state variables as a correctionterm obtained from an approximate solution of the real modal equationfor lower-order modes, and derive the dimension reduced model, and athird processing unit configured to calculate a frequency response byusing the dimension reduced model.
 7. A storage medium which has storedthe program according to claim
 6. 8. The vibration analysis methodaccording to claim 2, wherein, in the new type of complex modalanalysis, calculation of complex eigenvalues and complex eigenvectorsfor lower-order modes of linear state variables and realization of thecomplex eigenvectors are performed, and introduction of real modalcoordinates for physical coordinates of lower-order modes and derivationof real modal equations are performed.
 9. The vibration analysis methodaccording to claim 3, wherein, in the new type of complex modalanalysis, calculation of complex eigenvalues and complex eigenvectorsfor lower-order modes of linear state variables and realization of thecomplex eigenvectors are performed, and introduction of real modalcoordinates for physical coordinates of lower-order modes and derivationof real modal equations are performed.
 10. The vibration analysis methodaccording to claim 4, wherein, in the new type of complex modalanalysis, calculation of complex eigenvalues and complex eigenvectorsfor lower-order modes of linear state variables and realization of thecomplex eigenvectors are performed, and introduction of real modalcoordinates for physical coordinates of lower-order modes and derivationof real modal equations are performed.